In [2]:
%pylab inline
import random


Populating the interactive namespace from numpy and matplotlib

Metadynamics

Gaussian Potential Bias Analysis

The method of metadynamics seeks to explore a free energy landscape and cross high energy barriers using artificial bias potentials. This allows a system to pass through a highly unfavorable state without needing to quench the system. This is done by adding Gaussian potentials in deep energy wells, filling them, and allowing the system to "spill" into another metastable state.

$$ V_{bias} = G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$

In [3]:
figure(figsize=[6,2])
height = 5
center = 0
width = 2
x = linspace(-10,10,1000)
gaussian = height*exp(-(x-center)**2/(2*width)**2)
plot(x,gaussian)
title("Gaussian")


Out[3]:
<matplotlib.text.Text at 0x107fec110>

Adding a Gaussian potential to a free energy landscape

As the system moves forward in time, a Gaussian is added to the free energy surface $F(s)$ at the previous location.

$$ F'(r) = F(r) + \sum G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$

The plot below shows the effect of this bias potential on the surface. As more gaussians are added to the system, wells are filled and the MD simulation is able to sample the whole system.


In [3]:
figure(figsize=[10,5])
x = linspace(-1,1,1000)
height = .2
center = 0
width = .02
gaussian = height*exp(-(x-center)**2/(2*width)**2)
y = (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x)
y2 = (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x) + gaussian
fig = plt.figure(1)
plot(x,y, 'r')
plot(x,y2,'b')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')


Out[3]:
<matplotlib.text.Text at 0x1076ecd10>

Metadynamics simulation

This cell demostrates the behavior of an MD simulation that uses Metadynamics to sample all of surface. Notice there are a few key parameters that determine the computational expense and resolution of the results:

  1. Height of the Gaussian
  2. Width of the Gaussian
  3. Time step for each addition of the Gaussian

Depending on how rough the free energy landscape is, these parameters will need to be adjusted to smooth the surface and sample key metastable states.

In the example below, the red line describes the "real" free energy landscape (same arbitrary function for this example). The blue line is the sum of the bias potentials added to the free energy landscape. All the energy wells are sampled in this system and "filled" with artificial potentials.


In [4]:
figure(figsize=[10,6])

def y(x):
    return (3*x**10)+(2*x**9)-(x**8)+(3*x**7)-(2*x**6)-(5*x**5)-(6*x**4)+(x**3)+(7*x**2)-(x)
    
def add_gaussian(height, width, func):
    gaussian = height*exp(-(x-center)**2/(2*width)**2)
    return func(x) + gaussian

x = linspace(-1,1,1000)
surface = y(x)
y1 = surface
plot(x,surface, 'r')
height = .06
center = 0
width = .01
Gaussian = 0
for i in range(0,30000):
    center = random.uniform((-1*.2)+center,(1*.2)+center)
    gaussian = height*exp(-(x-center)**2/(2*width)**2)
    y2 = y1 + gaussian
    if max(y2) > 3:
        pass
    else:
        y1 += gaussian
        Gaussian += gaussian
plot(x,y1, 'b')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')


Out[4]:
<matplotlib.text.Text at 0x107907450>

Take the total artificial bias potential (sum of the Gaussians) and plot it's negative

$$ F(r) \approx - \sum G(s)e^{-\frac{(s-s_0)^2}{2\sigma ^2}} $$

Depending on how robust and resolved the metadynamics simulation is, the free energy landscape is approximately equal to this bias potential.


In [22]:
figure(figsize=[12,8])
plot(x,-Gaussian+3.1)
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')

figure()
figure(figsize=[12,8])
plot(x,-Gaussian+3.1)
plot(x,y(x), 'r')
title('Metadynamics Illustration')
xlabel('Reaction Coordinate')
ylabel('Free Energy')


Out[22]:
<matplotlib.text.Text at 0x10a4a2f90>
<matplotlib.figure.Figure at 0x10a499090>

Collective Variables

The importance of these type of free energy landscapes is highly dependent on the choice of collective variable's (a.k.a. reaction coordinates). There are three major criteria for CV's:

  1. They should clearly distringuish between the initial, final, and intermediate states (energy barrier is clearly seen).
  2. They should describe all the slow events that occur in the system
  3. They should have a limited number to increase speed of the simulation

The collective variables are the key to determining the mechanism for transitions between metastable states. Metadynamics takes a "course-grained" approach to sampling a system and finding these metastable states, focusing only on "interesting" variables that describe the dynamics of the system. These variables need to be chosen correctly to find results with meaning.